<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
                "http://www.w3.org/TR/REC-html40/loose.dtd">
<html>
<head>
  <title>Description of ADCP_DispCoef</title>
  <meta name="keywords" content="ADCP_DispCoef">
  <meta name="description" content="This program computes the longitudinal dispersion coefficient from ADCP">
  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
  <meta name="generator" content="m2html &copy; 2005 Guillaume Flandin">
  <meta name="robots" content="index, follow">
  <link type="text/css" rel="stylesheet" href="../../m2html.css">
  <script type="text/javascript">
    if (top.frames.length == 0) { top.location = "../../index.html"; };
  </script>
</head>
<body>
<a name="_top"></a>
<!-- ../menu.html VMT_4.0_dev --><!-- menu.html utils -->
<h1>ADCP_DispCoef
</h1>

<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="box"><strong>This program computes the longitudinal dispersion coefficient from ADCP</strong></div>

<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="box"><strong>function [k,kc,Eyc,Q] = ADCP_DispCoef(beddepth,travdist,vertdepth,downstvel,startDist,endDist,extrp,extend_to_banks,banktype); </strong></div>

<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="fragment"><pre class="comment">This program computes the longitudinal dispersion coefficient from ADCP
transects.</pre></div>

<!-- crossreference -->
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
This function calls:
<ul style="list-style-image:url(../../matlabicon.gif)">
<li><a href="fitLogLawV2.html" class="code" title="function [ustar,zo,cod] = fitLogLawV2(u,z,h)">fitLogLawV2</a>	This function fits a log law of the form u/u* = 1/kappa*ln(z/zo) to the given data</li><li><a href="nansum.html" class="code" title="function y = nansum(x,dim)">nansum</a>	NANSUM Sum, ignoring NaNs.</li></ul>
This function is called by:
<ul style="list-style-image:url(../../matlabicon.gif)">
</ul>
<!-- crossreference -->



<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="fragment"><pre>0001 <a name="_sub0" href="#_subfunctions" class="code">function [k,kc,Eyc,Q] = ADCP_DispCoef(beddepth,travdist,vertdepth,downstvel,startDist,endDist,extrp,extend_to_banks,banktype);</a>
0002 
0003 <span class="comment">%This program computes the longitudinal dispersion coefficient from ADCP</span>
0004 <span class="comment">%transects.</span>
0005 
0006 <span class="comment">%P.R. Jackson &amp; N.V. Reynolds, USGS, 11/16/10</span>
0007 
0008 <span class="comment">%Inputs:</span>
0009 <span class="comment">%</span>
0010 <span class="comment">%     beddepth  = Depth (in meters) from the water surface to the bed (n element vector)</span>
0011 <span class="comment">%     travdist  = Transverse distance (in m) across the river starting with (n element vector)</span>
0012 <span class="comment">%     vertdepth = Depths (in m) from the water surface to velocity bins in the water column (m element vector)</span>
0013 <span class="comment">%     downstvel = streamwise velocity distribution (+ for DS) in m/s (n x m array)</span>
0014 <span class="comment">%     startDist = distance (in m) from the first profile to the starting bank</span>
0015 <span class="comment">%     endDist   = distance (in m) from the last profile to the ending bank</span>
0016 <span class="comment">%     extrp     = (Binary) Set to 1 to extrapolate profiles to the bed and surface (log law); otherwise set to 0</span>
0017 <span class="comment">%     extend_to_banks  = (Binary) Set to 1 to extrapolate data to the banks; otherwise set to 0</span>
0018 <span class="comment">%     banktype  = (string) 'Triangular' or 'Square' (used in bank extrapolation)</span>
0019 
0020 <span class="comment">%Outputs:</span>
0021 <span class="comment">%     k  = Longitudinal dispersion coefficient (in m^2/s) (variable Ey)</span>
0022 <span class="comment">%     kc  = Longitudinal dispersion coefficient (in m^2/s) (Constant Ey)</span>
0023 <span class="comment">%     Ey = Transverse mixing coefficient (constant)</span>
0024 <span class="comment">%     Q  = Approximate discharge in m^3/s</span>
0025 
0026 
0027 
0028 <span class="comment">%% Basic parameters</span>
0029 nprof = length(beddepth);
0030 
0031 <span class="comment">%% Correct missing bed depths</span>
0032 <span class="comment">%figure(10); clf; plot(travdist,beddepth,'k-')</span>
0033 <span class="comment">% indx1 = find(isnan(beddepth));</span>
0034 <span class="comment">% indx2 = find(~isnan(beddepth));</span>
0035 <span class="comment">% if isnan(beddepth(1))</span>
0036 <span class="comment">%     beddepth(1) = beddepth(indx2(1))/2; %Fill end nans (not filled with interpolation)</span>
0037 <span class="comment">% end</span>
0038 <span class="comment">% if isnan(beddepth(end))</span>
0039 <span class="comment">%     beddepth(end) = beddepth(indx2(end))/2;  %Fill end nans (not filled with interpolation)</span>
0040 <span class="comment">% end</span>
0041 indx1 = find(isnan(beddepth));
0042 indx2 = find(~isnan(beddepth));
0043 beddepth(indx1) = interp1(travdist(indx2),beddepth(indx2),travdist(indx1));  <span class="comment">%fills the nans</span>
0044 
0045 <span class="comment">%figure(10); hold on; plot(travdist(indx1),beddepth(indx1),'r-')</span>
0046 
0047 <span class="comment">%% Fit each profile with a log law</span>
0048 
0049 <span class="comment">%[ustarXS,zoXS] = ADCP_ComputeXS_ShearVelocity(downstvel,vertdepth,beddepth)</span>
0050 <span class="comment">%</span>
0051 <span class="keyword">for</span> i = 1:nprof
0052     indx = find(~isnan(downstvel(i,:)));
0053     <span class="keyword">if</span> isempty(indx)
0054         <span class="keyword">if</span> i == 1 | i == nprof
0055             beep
0056             error(<span class="string">'Edge NaN Detected in Downstream Velocity'</span>)
0057         <span class="keyword">end</span>
0058         ustar(i) = nan;
0059         zo(i) = nan;
0060     <span class="keyword">else</span>
0061         z = beddepth(i)-vertdepth(indx);  <span class="comment">% Height above the bed</span>
0062         zmean(i) = nanmean(z);  <span class="comment">%Mean height above the bed</span>
0063         u = downstvel(i,indx);
0064         <span class="comment">%figure(10); clf; plot(downstvel(i,indx),beddepth(i)-vertdepth(indx))</span>
0065         <span class="comment">%[ustar(i),zo(i),ks,cod,upred,zpred,delta] = fitLogLaw(u',z,beddepth(i));</span>
0066         [ustar1(i),zo(i)] = <a href="fitLogLawV2.html" class="code" title="function [ustar,zo,cod] = fitLogLawV2(u,z,h)">fitLogLawV2</a>(u',z,beddepth(i));
0067     <span class="keyword">end</span>
0068 <span class="keyword">end</span>
0069 <span class="comment">%clean up problem values</span>
0070 ustar1 = real(ustar1);
0071 indx = find(isinf(ustar1)); <span class="comment">%turn if values to NAN</span>
0072 ustar1(indx) = nan;
0073 <span class="comment">%edge values</span>
0074 indx1 = find(isnan(ustar1));
0075 indx2 = find(~isnan(ustar1));
0076 <span class="keyword">if</span> isnan(ustar1(1))
0077     ustar1(1) = ustar1(indx2(1)); <span class="comment">%Fill end nans (not filled with interpolation)</span>
0078 <span class="keyword">end</span>
0079 <span class="keyword">if</span> isnan(ustar1(end))
0080     ustar1(end) = ustar1(indx2(end));  <span class="comment">%Fill end nans (not filled with interpolation)</span>
0081 <span class="keyword">end</span>
0082 indx1 = find(isnan(ustar1));
0083 indx2 = find(~isnan(ustar1));
0084 ustar1(indx1) = interp1(travdist(indx2),ustar1(indx2),travdist(indx1));  <span class="comment">%fills the nans</span>
0085 
0086 zo = real(zo);
0087 indx = find(isinf(zo)); <span class="comment">%turn if values to NAN</span>
0088 zo(indx) = nan;
0089 <span class="comment">%edge values</span>
0090 indx1 = find(isnan(zo));
0091 indx2 = find(~isnan(zo));
0092 <span class="keyword">if</span> isnan(zo(1))
0093     zo(1) = zo(indx2(1)); <span class="comment">%Fill end nans (not filled with interpolation)</span>
0094 <span class="keyword">end</span>
0095 <span class="keyword">if</span> isnan(zo(end))
0096     zo(end) = zo(indx2(end));  <span class="comment">%Fill end nans (not filled with interpolation)</span>
0097 <span class="keyword">end</span>
0098 indx1 = find(isnan(zo));
0099 indx2 = find(~isnan(zo));
0100 zo(indx1) = interp1(travdist(indx2),zo(indx2),travdist(indx1));  <span class="comment">%fills the nans</span>
0101 
0102 indx = find(zmean == 0); <span class="comment">%turn if values to NAN</span>
0103 zmean(indx) = nan;
0104 <span class="comment">%edge values</span>
0105 indx1 = find(isnan(zmean));
0106 indx2 = find(~isnan(zmean));
0107 <span class="keyword">if</span> isnan(zmean(1))
0108     zmean(1) = zmean(indx2(1)); <span class="comment">%Fill end nans (not filled with interpolation)</span>
0109 <span class="keyword">end</span>
0110 <span class="keyword">if</span> isnan(zmean(end))
0111     zmean(end) = zmean(indx2(end));  <span class="comment">%Fill end nans (not filled with interpolation)</span>
0112 <span class="keyword">end</span>
0113 indx1 = find(isnan(zmean));
0114 indx2 = find(~isnan(zmean));
0115 zmean(indx1) = interp1(travdist(indx2),zmean(indx2),travdist(indx1));  <span class="comment">%fills the nans</span>
0116 
0117 mzo = nanmedian(zo)  <span class="comment">%Median z0</span>
0118 <span class="comment">%figure(10); clf; plot(travdist,ustar','k',travdist,shearvel,'r',travdist,ustar2','b')</span>
0119 
0120 <span class="comment">%% Extrapolate top and bottom portions if required</span>
0121 <span class="comment">%figure(20); clf; contour(downstvel,'Fill','on','Linestyle','none');</span>
0122 test = downstvel;
0123 <span class="keyword">if</span> extrp  <span class="comment">%extrapolates to the bed and surface</span>
0124     zgridspc = nanmean(diff(vertdepth));
0125     topDist  = vertdepth(1)
0126     ntop   = floor(topDist/zgridspc);
0127     top_nodes   = linspace(0,topDist-zgridspc,ntop)';
0128     vertdepth   = [top_nodes; vertdepth];
0129     downstvel = [nan*ones(length(beddepth),ntop) downstvel];  <span class="comment">%preallocated space for new top velocities</span>
0130     <span class="comment">%figure(21); clf; contour(downstvel,'Fill','on','Linestyle','none');</span>
0131     
0132     <span class="keyword">for</span> i = 1:nprof
0133         indx = find(~isnan(downstvel(i,:)));
0134         
0135         <span class="comment">%figure(10); clf; plot(downstvel(i,:),vertdepth,'ko-'); set(gca,'YDir','reverse');</span>
0136         <span class="keyword">if</span> ~isempty(indx) <span class="comment">%all nans get no computations or extensions</span>
0137             
0138             top_nodes_flipped = beddepth(i) - top_nodes;
0139             u_top = ustar1(i)/0.41*log(top_nodes_flipped/zo(i));
0140             <span class="comment">%botDist = beddepth(i) - vertdepth(indx(end));</span>
0141             
0142             last_node = find(vertdepth &lt; beddepth(i));
0143             bot_nodes = vertdepth(indx(end)+1:last_node(end));
0144             
0145             <span class="comment">%nbot   = floor(botDist/zgridspc);</span>
0146             <span class="comment">%bot_nodes = linspace(vertdepth(indx(end))+zgridspc,vertdepth(indx(end))+botDist,nbot);</span>
0147 
0148             <span class="comment">%Flip the vertical coordinate and use log fit from above to extend</span>
0149             <span class="comment">%to boundaries</span>
0150 
0151             
0152             bot_nodes_flipped = beddepth(i) - bot_nodes;
0153             u_bot = ustar1(i)/0.41*log(bot_nodes_flipped/zo(i));
0154 
0155             <span class="comment">%u_bot(end) = 0.0;  %Log law is singular at z = 0, so reset to zero for no slip</span>
0156 
0157 
0158             <span class="comment">% Append to data</span>
0159             <span class="keyword">if</span> ~isempty(indx)
0160                 downstvel(i,indx(end)+1:last_node(end)) = u_bot;
0161             <span class="keyword">end</span>        
0162             <span class="comment">%downstvel(i,last_node(end)+1) = 0.0;  %Set to zero at bed</span>
0163             <span class="comment">%vertdepth(last_node(end)+1) = beddepth(i); %reset first invalid</span>
0164             <span class="comment">%bottom bin to bed depth</span>
0165             downstvel(i,1:ntop) = u_top;
0166         
0167         <span class="keyword">end</span>
0168 
0169         <span class="comment">%figure(10); hold on; plot(u_top,top_nodes,'r.',u_bot,bot_nodes,'b.'); set(gca,'YDir','reverse'); pause</span>
0170         
0171         clear indx last_node bot_nodes top_nodes_flipped bot_nodes_flipped u_top u_bot 
0172     <span class="keyword">end</span>
0173 <span class="keyword">end</span>
0174 <span class="comment">%test2 = downstvel(:,ntop+1:end) - test;</span>
0175 <span class="comment">%figure(22); clf; contour(test2,'Fill','on','Linestyle','none');</span>
0176 
0177 <span class="comment">%% Compute Depth-averaged velocities</span>
0178 
0179 <span class="comment">%This function computes the layer averaged mean of the downstream velocity over the measured depth range.</span>
0180 <span class="comment">%Assumes the data outside the depth range have been set to NaN.</span>
0181 <span class="comment">% P.R. Jackson, USGS 1-7-09</span>
0182 <span class="keyword">for</span> i = 1:nprof
0183     lam(i) = VMT_LayerAveMean(vertdepth,downstvel(i,:)');
0184 <span class="keyword">end</span>
0185 
0186 <span class="comment">%figure(10); clf; plot(travdist,lam,'r-') % Plot to check</span>
0187 
0188 <span class="comment">%lam = VMT_LayerAveMean(repmat(vertdepth,1,size(downstvel',2)),downstvel')';</span>
0189 lam = real(lam)';  <span class="comment">%remove any imaginary parts</span>
0190 indx = find(isinf(lam)); <span class="comment">%turn if values to NAN</span>
0191 lam(indx) = nan;
0192 indx1 = find(isnan(lam));
0193 indx2 = find(~isnan(lam));
0194 <span class="keyword">if</span> isnan(lam(1))
0195     lam(1) = lam(indx2(1))/2; <span class="comment">%Fill end nans (not filled with interpolation)</span>
0196 <span class="keyword">end</span>
0197 <span class="keyword">if</span> isnan(lam(end))
0198     lam(end) = lam(indx2(end))/2;  <span class="comment">%Fill end nans (not filled with interpolation)</span>
0199 <span class="keyword">end</span>
0200 indx1 = find(isnan(lam));
0201 indx2 = find(~isnan(lam));
0202 lam(indx1) = interp1(travdist(indx2),lam(indx2),travdist(indx1));  <span class="comment">%fills the nans</span>
0203 
0204 <span class="comment">%indx1 = find(isnan(lam))</span>
0205 <span class="comment">%figure(10); clf; plot(travdist,lam,'r-') % Plot to check</span>
0206 
0207 <span class="comment">% %Davide's method for LAM</span>
0208 <span class="comment">%</span>
0209 <span class="comment">% for i=1:sztravdist+1</span>
0210 <span class="comment">%     x(i)=nansum(downstvel(i,:));</span>
0211 <span class="comment">%</span>
0212 <span class="comment">% end</span>
0213 <span class="comment">%</span>
0214 <span class="comment">%</span>
0215 <span class="comment">% lamd=x'./beddepth;</span>
0216 <span class="comment">% lamd(1)=lamd(2)/2;</span>
0217 <span class="comment">% lamd(sztravdist+1)=lamd(sztravdist)/2;</span>
0218 <span class="comment">% %lamd=lamd';</span>
0219 <span class="comment">% figure(10); hold on; plot(travdist,lamd./3.281,'b-')</span>
0220 <span class="comment">% whos</span>
0221 <span class="comment">% return</span>
0222 
0223 <span class="comment">%% Recompute the shear velocity using the median z0 and LAM (lower noise)</span>
0224 <span class="comment">%shearvel = lam*0.41./log(zmean'./mzo);  %kappa = 0.41</span>
0225 <span class="keyword">if</span> 0  <span class="comment">%use the z0 from the normalized fit</span>
0226     shearvel = lam*0.41./log((beddepth/exp(1))./zoXS);  <span class="comment">%kappa = 0.41; %Use zmean = h/e following Sime etal 2007 as mean velocity occurs at h/e</span>
0227 <span class="keyword">else</span> <span class="comment">%uses the median of the individual profile fits</span>
0228     shearvel = lam*0.41./log((beddepth/exp(1))./mzo);  <span class="comment">%kappa = 0.41; %Use zmean = h/e following Sime etal 2007 as mean velocity occurs at h/e</span>
0229 <span class="keyword">end</span>
0230 <span class="comment">% replace zeros with nans</span>
0231 indx = find(shearvel == 0);
0232 shearvel(indx) = nan;
0233 <span class="comment">% Fill gaps</span>
0234 indx1 = find(isnan(shearvel));
0235 indx2 = find(~isnan(shearvel));
0236 shearvel(indx1) = interp1(travdist(indx2),shearvel(indx2),travdist(indx1));  <span class="comment">%fills the nans</span>
0237 
0238 <span class="comment">%figure(10); clf; plot(travdist,shearvel,'r')</span>
0239 travdist_orig = travdist;
0240 
0241 <span class="comment">%% Extrapolate to the banks (following WRII)</span>
0242 <span class="keyword">if</span> extend_to_banks
0243     ygridspc = nanmean(diff(travdist));
0244     nstart = floor(startDist/ygridspc)
0245     nend   = floor(endDist/ygridspc);
0246 <span class="comment">%Note: This fails for zero start or end distances</span>
0247     <span class="comment">%indx = find(~isnan(lam));</span>
0248 
0249     start_nodes = linspace(0,startDist-ygridspc,nstart);
0250     end_nodes   = linspace(travdist(end)+ygridspc,travdist(end)+endDist,nend) + startDist;
0251 
0252     <span class="keyword">switch</span> banktype
0253         <span class="keyword">case</span> <span class="string">'triangular'</span>
0254             startBed = interp1([start_nodes(1) startDist],[0 beddepth(1)],start_nodes);
0255             endBed   = interp1([(travdist(end)+startDist) end_nodes(end)],[beddepth(end) 0],end_nodes);
0256             startLAM = lam(1).*sqrt(startBed./beddepth(1));
0257             endLAM   = lam(end).*sqrt(endBed./beddepth(end));
0258             startSHV = interp1([start_nodes(1) startDist],[0 shearvel(1)],start_nodes);  <span class="comment">%Linear interpolation of shear velocities</span>
0259             endSHV   = interp1([(travdist(end)+startDist) end_nodes(end)],[shearvel(end) 0],end_nodes);
0260 
0261 
0262         <span class="keyword">case</span> <span class="string">'square'</span>
0263             errordlg(<span class="string">'Not implemented yet'</span>)
0264     <span class="keyword">end</span>
0265 
0266     travdist = [start_nodes'; travdist+startDist; end_nodes'];
0267     lam      = [startLAM'; lam; endLAM'];
0268     beddepth = [startBed'; beddepth; endBed'];
0269     shearvel = [startSHV'; shearvel; endSHV'];
0270     
0271     shearvel(1) = 0.1*shearvel(2);  <span class="comment">%Reset start and end u* to nonzero to keep INT2 from blowing up (due to zero Ey at edges)</span>
0272     shearvel(end) = 0.1*shearvel(end-1); 
0273 
0274     
0275     <span class="comment">%figure(10); hold on; plot(travdist,lam,'b-') % Plot to check</span>
0276     <span class="comment">%figure(10); clf; plot(travdist,beddepth,'b-') % Plot to check</span>
0277 
0278 <span class="keyword">end</span>
0279 
0280 <span class="comment">%% Compute the mean Cross-sectional area</span>
0281 XSarea = trapz(travdist,beddepth)
0282 
0283 
0284 <span class="comment">%% Compute bin areas</span>
0285 <span class="comment">%compute the midpoints</span>
0286 difftrav = diff(travdist);
0287 interim1 = 0.5*difftrav(1:end-1) + 0.5*difftrav(2:end);
0288 interim1 = [0.5*difftrav(1); interim1];
0289 bin_mp = cumsum(interim1);
0290 
0291 interim2 = [0; bin_mp; travdist(end)];
0292 delW = diff(interim2);
0293 
0294 <span class="comment">%figure(10); clf; plot(travdist,beddepth,'k.-');</span>
0295 
0296 <span class="comment">%Interpolate the bed depths at the midpoints</span>
0297 beddepth_mp = interp1(travdist,beddepth,bin_mp);
0298 
0299 <span class="comment">%figure(10); clf; plot(travdist,beddepth,'k.-',bin_mp,beddepth_mp,'r.'); %plot to check</span>
0300 
0301 <span class="comment">%Compute the areas</span>
0302 cumlarea = cumtrapz(interim2,[beddepth(1); beddepth_mp; beddepth(end)]);
0303 bin_areas = diff(cumlarea);
0304 <span class="comment">% Fill gaps</span>
0305 indx1 = find(isnan(bin_areas));
0306 indx2 = find(~isnan(bin_areas));
0307 bin_areas(indx1) = interp1(travdist(indx2),bin_areas(indx2),travdist(indx1));  <span class="comment">%fills the nans</span>
0308 
0309 
0310 <span class="comment">%figure(10); clf; plot(travdist,bin_areas,'k.-'); %plot to check</span>
0311 
0312 <span class="comment">%Check the computation</span>
0313 areaClosure = 100*(XSarea - sum(bin_areas))/XSarea
0314 
0315 
0316     <span class="comment">%Old area comp</span>
0317     <span class="comment">% sztravdist=size(travdist,1)-1;</span>
0318     <span class="comment">% szvertdepth=size(vertdepth,1);</span>
0319     <span class="comment">%</span>
0320     <span class="comment">% %midpoint to the left of y-component</span>
0321     <span class="comment">% midptly=(1:sztravdist+1);</span>
0322     <span class="comment">% midptly(1)=0;</span>
0323     <span class="comment">% midptly(2)=travdist(2)-(travdist(3)-travdist(2))/2;</span>
0324     <span class="comment">% midptly(3:sztravdist-0)=(travdist(3:sztravdist)+travdist(2:sztravdist-1))/2;</span>
0325     <span class="comment">% midptly(sztravdist+1)=travdist(sztravdist)-(travdist(sztravdist-1)-travdist(sztravdist))/2;</span>
0326     <span class="comment">%</span>
0327     <span class="comment">% %midpoint to the left of bed depth</span>
0328     <span class="comment">% midptld=(1:sztravdist+1);</span>
0329     <span class="comment">% midptld(1)=0;</span>
0330     <span class="comment">% midptld(2)=beddepth(2)-(beddepth(3)-beddepth(2))/2;</span>
0331     <span class="comment">% midptld(3:sztravdist-0)=(beddepth(3:sztravdist)+beddepth(2:sztravdist-1))/2;</span>
0332     <span class="comment">% midptld(sztravdist+1)=beddepth(sztravdist)-(beddepth(sztravdist-1)-beddepth(sztravdist))/2;</span>
0333     <span class="comment">%</span>
0334     <span class="comment">% %midpoint to the right of y-component</span>
0335     <span class="comment">% midptry=(1:sztravdist+1);</span>
0336     <span class="comment">% midptry(1)=travdist(2)-(travdist(3)-travdist(2))/2;</span>
0337     <span class="comment">% midptry(2:sztravdist-1)=(travdist(3:sztravdist)+travdist(2:sztravdist-1))/2;</span>
0338     <span class="comment">% midptry(sztravdist)=travdist(sztravdist)-(travdist(sztravdist-1)-travdist(sztravdist))/2;</span>
0339     <span class="comment">% midptry(sztravdist+1)=0;</span>
0340     <span class="comment">%</span>
0341     <span class="comment">% %midpoint to the right of bed depth</span>
0342     <span class="comment">% midptrd=(1:sztravdist+1);</span>
0343     <span class="comment">% midptrd(1)=beddepth(2)-(beddepth(3)-beddepth(2))/2;</span>
0344     <span class="comment">% midptrd(2:sztravdist-1)=(beddepth(3:sztravdist)+beddepth(2:sztravdist-1))/2;</span>
0345     <span class="comment">% midptrd(sztravdist)=beddepth(sztravdist)-(beddepth(sztravdist-1)-beddepth(sztravdist))/2;</span>
0346     <span class="comment">% midptrd(sztravdist+1)=0;</span>
0347     <span class="comment">%</span>
0348     <span class="comment">%</span>
0349     <span class="comment">% %area per bin</span>
0350     <span class="comment">% ameth2=(1:sztravdist+1);</span>
0351     <span class="comment">% ameth2(1)=(beddepth(1)+midptrd(1))*(midptry(1)-travdist(1))/2;</span>
0352     <span class="comment">% for i=2:sztravdist</span>
0353     <span class="comment">% ameth2(i)=(beddepth(i)+midptld(i)')*(travdist(i)-midptly(i)')/2+(beddepth(i)+midptrd(i)').*(midptry(i)'-travdist(i))./2;</span>
0354     <span class="comment">% end</span>
0355     <span class="comment">% ameth2(sztravdist+1)=((beddepth(sztravdist+1)+midptld(sztravdist+1))*(travdist(sztravdist+1)-midptly(sztravdist+1)))/2;</span>
0356 
0357 <span class="comment">%figure(10); clf; plot(travdist,bin_areas,'k.-',travdist,ameth2,'r.-')</span>
0358 
0359 <span class="comment">%Compute the depth for each subarea</span>
0360 d = bin_areas./delW;
0361 <span class="comment">% Fill gaps</span>
0362 indx1 = find(isnan(d));
0363 indx2 = find(~isnan(d));
0364 d(indx1) = interp1(travdist(indx2),d(indx2),travdist(indx1));  <span class="comment">%fills the nans</span>
0365     <span class="comment">% %Old depth comp</span>
0366     <span class="comment">% d2=(1:sztravdist+1);</span>
0367     <span class="comment">% d2(1)=ameth2(1)/midptry(1);</span>
0368     <span class="comment">% d2=ameth2./(midptry-midptly);</span>
0369     <span class="comment">% d2(sztravdist+1)=ameth2(sztravdist+1)./(-midptly(sztravdist+1)+travdist(sztravdist+1));</span>
0370     <span class="comment">% d2 = d2';</span>
0371     
0372 <span class="comment">%figure(10); clf; plot(travdist,d,'k.-')</span>
0373 
0374 
0375 <span class="comment">%% Compute the Approximate Discharge</span>
0376 Q = <a href="nansum.html" class="code" title="function y = nansum(x,dim)">nansum</a>(lam.*bin_areas)
0377 <span class="keyword">if</span> Q &lt; 0
0378     Q = -Q;
0379     lam = -lam;
0380     shearvel = -shearvel;
0381 <span class="keyword">end</span>
0382 
0383 <span class="comment">%% Compute the mean cross-sectional velocity</span>
0384 Ua = Q./XSarea
0385 
0386 <span class="comment">%% Compute the LAM velocity deviation</span>
0387 uprime = lam - Ua;
0388 <span class="comment">%figure(10); clf; plot(travdist,uprime.*d,'b-')</span>
0389 
0390 <span class="comment">%% Compute the transverse mixing coefficients (variable)</span>
0391 B = travdist(end) - travdist(1)  <span class="comment">%Width</span>
0392 <span class="keyword">if</span> 1
0393     Ey = 0.21*shearvel.*d;  <span class="comment">%Elder approach</span>
0394     <span class="comment">% Fill gaps</span>
0395     indx1 = find(isnan(Ey));
0396     indx2 = find(~isnan(Ey));
0397     Ey(indx1) = interp1(travdist(indx2),Ey(indx2),travdist(indx1));  <span class="comment">%fills the nans</span>
0398 <span class="keyword">else</span>
0399     <span class="comment">% Compute the value following Shen et al. 2010</span>
0400     Ustar = trapz(travdist,shearvel)/B  <span class="comment">%Cross-sectionally averaged shear velocity</span>
0401     H = XSarea/B
0402 
0403     theta = 0.145 + (1/3520)*(Ua/Ustar)*(B/H)^1.38
0404 
0405     Ey = theta*Ustar*d;
0406 
0407     <span class="comment">% Fill gaps</span>
0408     indx1 = find(isnan(Ey));
0409     indx2 = find(~isnan(Ey));
0410     Ey(indx1) = interp1(travdist(indx2),Ey(indx2),travdist(indx1));  <span class="comment">%fills the nans</span>
0411 <span class="keyword">end</span>
0412 
0413 Eyc = trapz(travdist,Ey)/B;  <span class="comment">%Constant value of tranverse mixing coefficient</span>
0414 <span class="comment">%figure(10); clf; plot(travdist,Ey,'b-')</span>
0415 
0416 <span class="comment">%% Compute the integrals (variable Ey)</span>
0417 <span class="comment">% Integral 1</span>
0418 <span class="comment">%indx = find(~isnan(uprime));</span>
0419 INT1 = cumtrapz(travdist,uprime.*d);
0420 <span class="comment">%INT1 = cumtrapz(travdist(:),uprime(:).*d(:));</span>
0421 
0422 <span class="comment">% Integral 2</span>
0423 INT2 = cumtrapz(travdist,INT1./(Ey.*d));
0424 <span class="comment">%INT2 = cumtrapz(travdist(:),INT1./(Ey(:).*d(:)));</span>
0425 
0426 <span class="comment">% Integral 3</span>
0427 INT3 = cumtrapz(travdist,INT2.*uprime.*d);
0428 <span class="comment">%INT3 = cumtrapz(travdist(indx),INT2.*uprime(indx).*d(indx));</span>
0429 
0430 <span class="keyword">if</span> 0
0431     figure(10); clf;
0432     title(<span class="string">'Variable Ey'</span>)
0433     subplot(3,1,1)
0434     plot(travdist,INT1,<span class="string">'b-'</span>)
0435     ylabel(<span class="string">'INT1'</span>)
0436     subplot(3,1,2); plot(travdist,INT2,<span class="string">'b-'</span>)
0437     ylabel(<span class="string">'INT2'</span>)
0438     subplot(3,1,3); plot(travdist,INT3,<span class="string">'b-'</span>)
0439     ylabel(<span class="string">'INT3'</span>)
0440 <span class="keyword">end</span>
0441 
0442 <span class="comment">%% Compute and return the longitudinal dispersion coefficient</span>
0443 
0444 k = -INT3(end)/XSarea;
0445 <span class="comment">%percentdiff=100*(2111.7-k)/(2111.7);</span>
0446 disp(<span class="string">' '</span>)
0447 disp([<span class="string">'The longitudinal dispersion coefficient(K)is '</span> num2str(k) <span class="string">' m^2/s'</span>])
0448 disp(<span class="string">'          with a VARIABLE transverse mixing coefficient(Ey) '</span>)
0449 <span class="comment">%disp(' ')</span>
0450 <span class="comment">%disp(['Davide ended up with 2111.7 m^2/s. Which is a ' num2str(percentdiff) '% difference.  '])</span>
0451 
0452 <span class="comment">%% Compute the integrals (constant Ey)</span>
0453 
0454 <span class="comment">% Integral 2</span>
0455 INT2c = cumtrapz(travdist,INT1./(Eyc.*d));
0456 <span class="comment">%INT2 = cumtrapz(travdist(:),INT1./(Ey(:).*d(:)));</span>
0457 
0458 <span class="comment">% Integral 3</span>
0459 INT3c = cumtrapz(travdist,INT2c.*uprime.*d);
0460 <span class="comment">%INT3 = cumtrapz(travdist(indx),INT2.*uprime(indx).*d(indx));</span>
0461 
0462 <span class="keyword">if</span> 0
0463     figure(11); clf; title(<span class="string">'Constant Ey'</span>)
0464     subplot(2,1,1); plot(travdist,INT2c,<span class="string">'b-'</span>)
0465     ylabel(<span class="string">'INT2'</span>)
0466     subplot(2,1,2); plot(travdist,INT3c,<span class="string">'b-'</span>)
0467     ylabel(<span class="string">'INT3'</span>)
0468 <span class="keyword">end</span>
0469 
0470 <span class="comment">%% Compute and return the longitudinal dispersion coefficient</span>
0471 
0472 kc = -INT3c(end)/XSarea;
0473 <span class="comment">%percentdiff=100*(2111.7-kc)/(2111.7);</span>
0474 disp(<span class="string">' '</span>)
0475 disp([<span class="string">'The longitudinal dispersion coefficient(K)is '</span> num2str(kc) <span class="string">' m^2/s'</span>])
0476 disp(<span class="string">'          with a CONSTANT transverse mixing coefficient(Ey) '</span>)
0477 
0478 <span class="comment">%% Plot a combined plot</span>
0479 figure(8); clf
0480 subplot(3,1,1); plot(travdist,beddepth,<span class="string">'k-'</span>,travdist,d,<span class="string">'r-'</span>); ylabel(<span class="string">'Bed Depth (m)'</span>)
0481 <span class="comment">%subplot(4,1,2); plot(travdist_orig,zo,'k-'); ylabel('z0 (m)')</span>
0482 subplot(3,1,2); plot(travdist,lam,<span class="string">'k-'</span>); ylabel(<span class="string">'Depth-Averaged Vel. (m/s)'</span>)
0483 subplot(3,1,3); plot(travdist,shearvel,<span class="string">'k-'</span>); ylabel(<span class="string">'u* (m/s)'</span>)
0484 xlabel(<span class="string">'Transverse Distance (m)'</span>)
0485 figure(9); clf
0486 subplot(3,1,1); plot(travdist,bin_areas,<span class="string">'k-'</span>); ylabel(<span class="string">'Bin Areas (m^2)'</span>)
0487 subplot(3,1,2); plot(travdist,uprime,<span class="string">'k-'</span>); ylabel(<span class="string">'uprime (m/s)'</span>)
0488 subplot(3,1,3); plot(travdist,Ey,<span class="string">'k-'</span>); ylabel(<span class="string">'Ey (m^2/s)'</span>)
0489 xlabel(<span class="string">'Transverse Distance (m)'</span>)
0490 
0491 
0492 
0493 
0494 
0495 
0496 
0497</pre></div>
<hr><address>Generated on Thu 21-Mar-2013 09:32:01 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" target="_parent">m2html</a></strong> &copy; 2005</address>
</body>
</html>